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Abstract 

In order to simulate the mechanical behavior of large structures assembled from thin composite panels, we 
propose a coupling technique which substitutes local 3D models for the global plate model in the critical zones 
where plate modeling is inadequate. The transition from 3D to 2D is based on stress and displacement distributions 
associated with Saint-Venant problems which are precalculated automatically for a simple 3D cell. The hybrid 
plate/3D model is obtained after convergence of a series of iterations between a global plate model of the structure 
and localized 3D models of the critical zones. This technique is nonintrusive because the global calculations can be 
carried out using commercial software. Evaluation tests show that convergence is fast and that the resulting hybrid 
model is very close to a full 3D model. 

keywords: Plate-3D coupling; mixed dimensionality; nonintrusive coupling; laminated composites. 


1 Introduction 

This paper deals with nonintrusive techniques for the simulation of structures assembled from large flat composite 
panels. In many such structures, one dimension is much smaller than the other two, so they can be modeled as plates. 
However, plate theories alone cannot take into account the 3D stress concentration states which occur near edges, 
geometric discontinuities or zones subjected to concentrated loads. These phenomena are important in the case of 
laminated composites because they are often responsible for the initiation and propagation of damage, e.g. in the case 
of delamination [20]. Therefore, a natural approach to dealing with large composite structures is to attempt, whenever 
necessary, to couple plate or shell models with 3D models [18]. 

The main topic of this paper is the development and evaluation of a type a coupling between a plate model and 
a 3D model which could be applied directly within the nonintrusive framework proposed in [23, 22] as well as in [5] 
for dynamics. Essentially, the nonintrusive coupling technique consists, starting from an initial global calculation, in 
alternating between local calculations with prescribed displacements at the boundaries and global calculations which 
include corrective loads in order to reduce the imbalance between the models. Because the nonintrusive coupling 
technique used in the paper was designed for models with the same dimensionality, the main problem addressed here 
is the recovery of the 3D quantities associated with the 2D model. 

Overall, the question of the link between a plate model and its 3D counterpart is of such practical and theoretical 
importance that it has given rise to many publications. This introduction focuses more specifically on the references 
of interest in the literature that we tried to incorporate into our formulation. 

The quality of a plate model is usually assessed by comparison with a corresponding 3D reference model. This 
has led to a first group of works focusing on the development of a posteriori error estimators, initially neglecting edge 
effects [28, 14, 47, 29, 30], then taking them into account [31]. Thus, by the end of the 1980’s, the performance of the 
classical Kirchhoff-Love and Reissner-Mindlin plate and shell theories in 3D elasticity was relatively well-understood. 
Another group of works concerned the technique of asymptotic developments with the thickness as the small parameter 
[19, 24, 12, 17, 16, 6, 11]. In addition, using the same approach as Reissner [40, 42] in statics and Mindlin in dynamics 
[41], a shear factor was often introduced [49] in order to correct what would have been deduced from a pure kinematic 
assumption by taking into account some a priori information on the three-dimensional stress state. This mixed nature 
of the link between plate theory and 3D theory is the reason why the construction of improved theories [35, 48, 34, 39, 8] 
and associated finite elements [9, 2] relies heavily on mixed formulations. 

Regardless of the method, the assessment of a plate model involves the comparison between a 3D “Saint-Venant” 
or long-wavelength solution [32, 25] and the solution obtained from plate theory. A key aspect of the comparison. 
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especially in the case of composites, concerns the distribution of the stress field [38, 15]. Therefore, a basic technique 
which is used, sometimes implicitly, when estimating the quality of plate theories is the recovery of equilibrated 3D 
solutions by integrating the equilibrium equations through the thickness [43, Ij. The application of this technique to 
finite element calculations requires the recovery of 2D fields which satisfy the 2D equilibrium equations exactly. The 
use of such a stress recovery technique [37, 13] makes the implementation within a commercial program difficult. A 
comparative analysis of various models and techniques to evaluate a priori or a posteriori out-of-plane stresses was 
presented in [7]. An improved hybrid post-processing procedure was proposed in [10] along with a new displacement- 
based element. This procedure, which works element-wise and, thus, is easier to implement, was extended to first-order 
theories in [27] to recover both displacements and stresses. 

In order for the 3D approximation derived from plate theory to be valid everywhere, it should be completed by 
adding an approximation of the edge effects. A key issue is that the problem deduced from the plate theory residual 
at the edges must lead to a localized problem [31, 33, 1]. Otherwise, even using what one calls refined theories, the 
plate’s interior solution would be meaningless for the orders of magnitude considered. This may be one of the reasons 
why the Reissner-Mindlin theory is the dominant approach today, even for composite structures. Another reason is, 
of course, that because it involves only first-order derivatives it is much easier to build finite elements based on this 
theory than, for example, on Kirchhoff-Love’s. However, the Reissner-Mindlin theory includes an artificial edge effect 
which, if not handled properly, could pollute the recovery of the 3D solution near the edges [44]. 

In this context, our approach to nonintrusive coupling must comply with several constraints. The first constraint 
is that our method involves two overlapping models, a 2D (plate or shell) model which is defined over the whole 
domain to be analyzed and remains unchanged throughout the analysis, and a refined 3D model in which structural 
details and complex nonlinearities such as contact or delamination can be taken into account. The second constraint 
is that we need to use results obtained from commercial finite element codes which, usually, are based on the Reissner- 
Mindlin model and do not provide recovery techniques. The implementation of such recovery techniques would require 
rather complex manipulations with simultaneous access to data from several elements, which would annihilate the 
very purpose of nonintrusive techniques. Moreover, the plate model is usually less than optimal compared to its 3D 
counterpart in that the plate quantities which would be deduced from the calculation of the full 3D model would be 
somewhat different from those obtained using the plate model, even far away from the edges. That is the reason why 
[37, 51] and others have used an asymptotic formulation up to the second order to define a quasi-optimal version of 
Reissner-Mindlin’s composite plate theory by means of an optimization procedure. Then, based on these accurate 
plate quantities, a recovery technique can be used to obtain the 3D displacement, strain and stress fields as functions 
of 2D variables. 

The paper is organized as follows: Section 2 shows how relevant 3D quantities can be generated for any stacking 
sequence in order to recover 3D stress and displacement bases. Then, in Section 3, these bases are used to define 
hybrid models. In Section 4, a first analysis of the results enables us to identify the most favorable coupling technique, 
which is tested on a more representative example in Section 5. 


2 Numerical recovery of the Saint-Venant stress and warping 

In this section we investigate the association of 3D quantities (traction and warping) with composite plate calculations. 


2.1 Notations and equations 

Let us consider a classical 3D problem defined in domain Let / be the volume force, g the traction prescribed 
over and Ud the prescribed displacement over the complementary part The unit vectors of an orthogonal 

spatial basis are denoted by (e^, ey, Cz). The linear elastic behavior is given by Hooke’s tensor H. The symmetric 
part of the displacement gradient and the Cauchy stress tensor are designated respectively by e(u) and a. 

We introduce the affine space of the admissible displacements W(D): 

W(D) = {w G i7^(D), u = Ujj on 


whose associated vector space is Uo{fl). 

The 3D problem is: 

Find u G W(D) such that, Vm* G Wo(D), 


[ 2L • ^{yt)dx = / / ■ ytdx + f g ■ u*dS 

Jq Jq J 

CT = H : e{u) 


( 1 ) 


Regarding the plate model, we have D = a; x [—h, h]. 

Using the notation ~ to designate in-plane quantities (x and y coordinates), the first assumption concerns the shape 
of the plate’s displacement: 

v=y + VzOz (2) 
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where v and Vz are respectively the in-plane and out-of-plane displacements of the mid-surface and 9 denotes the 
rotation of the section, all of which are functions defined in lo. Operator A is the cross product defined by: 


fxi\ /yi\ /x2y3-x3y2\ 

2^2 A h/2 = Xsyi - Xiy3 

yxs/ yys/ \x1y2-x2y1J 


This leads to the introduction of the classical plate strains: 


dun = 


7® = 

— \sym 


—2 

sym Vz,z 

'^X d~'^y , 3 


with 


I = xiv) + zgii) 

^i^z) + ^ A 






' n ^y,y 

_4. 


Using notation () for the integration through the thickness, the plate’s stress tensors are: 

IV = (ct) M = (za) Q = (a^) 


(3) 


(4) 


The second assumption of the plate model is that the peeling stress cr 2 z is negligible compared to the other 
components of the stress tensor. This assumption simplifies the expression of the plate’s behavior as a function of the 
plane stress Hps and out-of-plane stress B. (In 3D, = Be^.) 

For a plate in equilibrium, the constitutive relations are: 

K = (tp«) : giv) 

K = : gii) (5) 

Q = k{B) : e^ivz,9) 

where k is a shear correction factor. (In this paper, we use the correction factor from [3].) 

The admissible plate displacement field is: 

U^{ui) = ^iv,9,Vz) G 7V^(a;)®,which satisfies Dirichlet’s BCs| (6) 

and the plate problem becomes^: 

Find (v, 9 ,Vz) G {uj) such that V(v *,0 ,v*) G (uj) 

< [ (R-l{r)+K-xii*) + Q-U<l*))dx= [ l-v*dx+ [ g-v*dS (7) 

JUJ ^ ' JQ J 

^K={%s):±{v), :|(0), Q = k(B) :e,(n2,0) 


Remark 1. The Reissner-Mindlin kinematic assumption was chosen because it is the most commonly used in com¬ 
mercial software. This model provides an estimation of the transverse .shear forces, but, unfortunately, it leads to 
perturbations of the generalized forces near the edges [44 ]■ 

Remark 2 (Notations). From now on, we will use the following notations: 
for the kinematics: 

plate DOFs V = [vx Vy 
plate kinematics ® = (e^x 
and for the generalized forces: 

B = {Nxx Nxy 
SO that Fi ;5 = (ct : (^1.5 
^6:8 = : @2,4,5 

Although these notations seem biased toward the (e^,,^^) coordinate system (especially toward direction Cx since the 
components related to this normal direction are listed first), the following calculations restore the normal-independent 
property of the method. 

In order to deal with interfaces which are not aligned with the axes, we introduce the following notations: for 
an interface whose local coordinate system is {n,t,e^) (n being the normal vector and t the tangent vector), de¬ 
notes the angle such that cos{(j)) = n • Cx — The kinematic basis of the plate in the local coordinate sys¬ 

tem is = {n,t,—zn,zt,efi). The associated DOFs are = {v„,Vt,9n,9t,Vz) and the stress components are 

= [Nnn Nnt Mnn Qn Ntt Mtt Qt) ■ 

^for the sake of simplicity, we did not develop the right-hand-side in terms of plate components 


Ox Oy Vz)' 


z^x e,) 


SO 


that V = E 




AIxx AIxy Qx Tlyy Alyy Qy) 
'^x)) 

®ey)) 


( 8 ) 
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2.2 Recovery of the Saint-Venant traction fields 

Since plate solutions are relevant only for problems with long characteristic lengths of variation, we try to associate 
one 3D Saint-Venant stress held (r.) with each plate stress component. 

In order to do that, we consider the simplest 3D problems whose solutions are Saint-Venant helds, which consist 
in a sufficiently large square plate with the same material, the same stacking sequence and the same hnite elements 
as the local model, subjected to 8 components of generalized forces (3 in tension, 3 in bending and 2 in shearing). 
This leads to 8 basic problems with well-chosen regular loads, each activating one generalized force as dehned in 
Table 1 and illustrated in Figures 1-5. We call these problems “cell problems” because of the analogy of the method 
with micro-macro approaches. Here, we are less interested in the average response than in the distribution of the 3D 
mechanical quantities through the thickness in the center of the plate, which represents a numerical “3D recovery” of 
the plate quantities. 


j 

loading 

main load 

figure 

1 

<^xx = “1 over clDo, a^x = 1 over OFIl 

Nxx 

1 

2 

cTxy = ~1 overdflL <^xy = 1 overclDo 
cTxy = ~1 over9Da fjxy = I over9r2 a 

Nxy 

2 

3 

cTxx = z over 9f2o, a^x = z over BVIq 

Mxx 

3 

4 

cTxy = —z overSflo axy = z over d fir 
o'xy = z over dFl-a cfxy = — 2 ; over9Da 

Mxy 

4 

5 

CTrcz = 1 overcAIi 

Qx 

5 

6 

same as Problem 1, but rotated 90° 

^yy 


7 

same as Problem 3, but rotated 90° 

Myy 


8 

same as Problem 5, but rotated 90° 

Qy 



Table 1: Construction of the Saint-Venant solutions 


From each calculation j € |1,8|, one can extract the distribution of the 3D stress tensor (a^) through the thickness 
in the center of the plate, i.e. sufficiently far away from any edge effect. This is illustrated in Figure 6 after 
FE discretization. From this stress tensor, using Equation (8), one can post-process the generalized forces (F^) of 
the plate. These test cases were designed so that Fj is dominant in the experiment (although coupling among 
generalized forces exists naturally, if only because of the conservation of moments which involves shear forces). 

Because of the linearity of the problem, any Saint-Venant stress state is a combination of the extracted (a.). The 
(t^) are constructed such that the contributions of the generalized forces are uncoupled: 


> FW.=cr. 

J=i =J 

i=l 






(9) 


where brackets denote matrices. Then, the (r .) are obtained by inversion of the 8x8 matrix [F*]. 

In other words, the (t .) belong to the space of the Saint-Venant stresses and satisfy the orthogonality relation: 


1 ^ ^ p y for n = a: and i = fc e |I, 51 

—i '^—^ or for u = 2 /and (i, fc) e {(6,1), (7,3), (8, 5)} (10) 

0 = (Zj : ® e„)) in all other cases 

Figure 7 shows typical distributions of the nodal forces through the thickness, both for an isotropic plate and for 
4 orthotropic plies. (A complete description of the material will be presented in Section 4.) In the isotropic case, we 
recover the analytical functions of [36]. 

Remark 3. The (r^) were designed for loads which are aligned with axes {e^,ey). If necessary, one can substitute 
a tensor family {rf) in a rotated direction: given a local coordinate system {n,t), let n = n^e,^ + '^y&.y ■so that 
t = —nyC,,. + UxCy. Using linearity, the Saint-Venant solution for a load which is aligned with the local frame is: 



= -n,,nyZ,^ +inl- nl)z^^ + UyU^^z 


=N„ 


= nlZf^ - n,,nyz + 


= + ^vL. 


■Qy ’ =Q 


XT, = -nyZ ny,z. 


( 11 ) 


Obviously, the same transformation can be used for moments as well as for tensions. 
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Figure 6: Extraction of the stress profiles and nodal forces 


Figure 5: Cell problem, loading Q 
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Remark 4. The quality of the recovered Saint- Venant stress fields can be improved iteratively: the calculated stress 
fields are applied as traction boundary conditions to the cell problem, resulting in a problem with even more localized 
edge effects and a better stress distribution through the thickness in the center of the domain. Figure 8 shows that the 
stress field converges very rapidly because each correction is an order of magnitude smaller than the previous one. 




(a) Shear loading, xz component 


(b) Shear loading, zz component 


Figure 7: Distribution of the nodal forces through the thickness for isotropic and orthotropic materials 



Figure 8 : The norm of the correction to the Saint-Venant stress fields after each iteration 


2.3 Recovery of the Saint-Venant warping vectors 

Now let us try to define relevant information concerning the displacement field. Since the plate’s kinematics (rigid 
section) is unrealistic and would lead to significant errors [ 1 ], we propose to correct this kinematics using warping 
vectors {i.e. 3D displacements which are allowed to vary through the thickness) deduced from the Saint-Venant 
auxiliary problems. 

From the eight numerical experiments, we obtain the stress states {a.) and the displacement profiles (uN. By 

inverting the matrix of the generalized forces of the plate [F*], one obtains the family (t.) of the uncoupled stresses. 
The same matrix can be used to define the displacements C/j associated with pure Saint-Venant loads: 

[C/i ... lU] = [ui (12) 

Each of these displacements consists of a plate part and a warping part. In order to separate the two contributions, 
we assume that the plate displacements produce the same work as the 3D displacements in the Saint-Venant tractions. 
At each point of interface 7 , one has: 

V? = (Z, :(C/.®eJ) (*Gll,8l, fee [1,51) 

^ (13) 

k=l 

Figure 9 presents some of the warping functions extracted from the numerical 3D recovery. The result is in the 
form of classical zigzags [ 8 ] with problem-specific shapes. 

Remark 5. The extraction of the five kinematic plate components is carried out using only the first five Saint-Venant 
stress tensors e,,.: (r^ • ex)k=i :5 the normal direction. Nevertheless, this definition is sufficient for the warping to 
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be well-defined (apart from negligible terms) and orthogonal to all eight Saint-Venant stress tensors for any normal 
direction. Indeed, from [1 ], we know that for Saint- Venant problems the plate displacement corrected by warping in 
direction results in 3D stresses which are 0(h/L)-accurate, 2h being the thickness of the plate; additional warping 
in the plane leads to 0(h‘^/L"^)-accuracy. This leads to the following orders of magnitude (which, in practice, turn out 
to be very pessimistic because the numerical measurements are much smaller): 


((Zj • - Z 2 • Gy) • Uf) = ((^2 • - Zg • ey) ■ Uf) 

((Z 3 • - Z 4 • Sy)) = - Z 7 • Sy) ■ Uf) 

((Zg • - Zg • ey) • Uf) 


o{h^lL^){Kr.l) ' 

o{hyL^)iM,--g) > 

o{h/L){Q^-i„f) _ 


Vze [1,81 
VjG {1,2,6} 
Vfc G {3,4,7} 
V/ G {5,8} 


(14) 


These relations imply that the warping vectors are independent of the axes chosen and produce zero work in all Saint- 
Venant stress fields, even those that are rotated (jf^) because of linearity. 




(a) Shear warping, normal component 


(b) Shear warping, tangential component 


Figure 9: The extracted warping functions for isotropic and orthotropic materials 


3 Hybrid models 

In this section, we focus on the definition of hybrid models for thin structures using plate theories wherever appropriate, 
but keeping 3D models when necessary (see Figure 10). In our context of using commercial software, only less-than- 
optimal plate theories are available and, thus, a direct plate/3D connection would be error-prone. Besides, the analysis 
of an industrial structure often begins with a global plate model, so we want to take advantage of this model to the 
greatest possible extent and supplement it with local 3D modeling only where necessary. 

By convention, we use lowercase for the plate domain and for the interfaces (w, 7 ) and uppercase for their 3D 
counterparts (D = w x \—h, h], F = 7 x [—h, h]). We distinguish three zones: the zone of interest I, the buffer zone B 
and the complementary zone C (which exists only in the plate model). 

Then, the hybrid 3D plate model is defined as the limit of an iterative process, illustrated in Figure 11, which 
alternates between global plate calculations (domain w = w/ Uws Uwc) and local 3D calculations (domain D/ U D^). 






Figure 10: The reference 3D model and the hybrid model 


7 




















corrective 



First, a global plate calculation is used to define the boundary conditions for the local 3D problems along 7 c (the 
zooming step). As explained in Section 3.1, the plate-to-3D recovery uses the stress and warping distributions through 
the thickness deduced from the cell problems of Section 2. 

After the local calculation, the imbalance, in a plate sense, between the plate model and the 3D model is evaluated 
along the inner line 7 / (Section 3.2). Then, if necessary, a global plate calculation under a corrective loading applied 
along 7 / is carried out (Section 3.3). 

Beside the dimensionality transition problem, our approach differs from the algorithm of [23] in that we propose 
to distinguish interface 7 c (which serves to prescribe the boundary conditions of the zone of interest) from interface 
7 / (which is used to apply the global correction) and introduce a buffer zone Ds in-between. For the descent steps 
described in Sections 3.1.3 and 3.1.2, it is possible to choose 7 c = 7 /, but spurious edge effects may occur near Fc. 
The buffer zone Ds is used to let the evanescent components activated by the boundary conditions become negligible 
before entering the zone of interest D/ itself. In a sense, the plate model and the 3D model are made equivalent in 
the buffer zone, which makes this method different from other approaches based on a volume connection zone, such 
as the transition element method [21, 45] or the Arlequin method [4]. Besides, the implementation of the method is 
simpler thanks to its compatibility with nonintrusive approaches. 


3.1 The zooming step 

In this section, we assume that a plate problem has been solved and that all its kinematic and static quantities are 
known along line 7 c. Let (]Fi)ie|i,8] and (Vfc)feg|i^ 5 j denote respectively the 8 static plate components and the 5 
displacement components at the interface, and let s be the curvilinear abscissa of interface 70 . Thus, the position of 
a point of Fq is given by (s, z) and the classical 3D plate displacement is v^{s, z) = ^fe(-^)V/c(s)- 

The objective of the zooming step is to define relevant boundary conditions along Fq = 7 q x [—h,h] in order to 
solve the 3D problem in the extended zone of interest ilj U Ds (see Figure 10). In order to do that, we use the stress 
and warping distributions (r^) and (w^)s calculated for the cell problems of Section 2. 


3.1.1 Pure traction descent 

At each point of interface Fq, the following traction distribution is prescribed: 
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g{s,z) -nis) = 

2=1 


(15) 


Then, the 3D problem becomes: 

Find u e U{Vti) such that My* G Uq{VLi) 


I (H : |(u)) 


:e{u*)dx = J f ■ u*dx + J g ■ u*dS + J Fi(s)(r^.(z) ■ njs)) ■ u*dsdz 


(16) 


Of course, this technique works only for domains of interest possessing sufficient Dirichlet conditions. Therefore, it 
will be evaluated only briefly. 


3.1.2 Lagrangian descent 

The Lagrangian technique is a kinematic coupling method based on an energy equivalence which is similar to that 
proposed in [36, 46] for isotropic plates with known analytic solutions in order to calculate the distribution of Saint- 







Venant stresses (zf)- Here, since we are considering stacks of orthotropic plies, no such formula is available and, 
therefore, we propose to use the numerically calculated (zf)- 

The zooming step consists in prescribing the 3D displacement at the interface which develops the same amount of 
work as the given plate displacement with the Saint-Venant stresses (zf)- 

Find u € U(Qi), (Afc) € A( 7 c)® such that Vu* G 

f H : e{u) : €^{u*)dx = ( f ■ u*dx + f g ■ u*dS 
«/ 0 / Jq,j 

5 

+ f Afe(s) f u* ■ zf(z) ■ n(s) dz ds 
+ X! / ^f(s}(zf(z)-n(s))-u*dzds 

i=6 

[ \Us) [{u-v^)-zt{z)-]l{s)dzds = 0, Vfc e |l,5l,VAfc G A( 7 c) 

Jyc J z 

where A( 7 c’) is the space of the Lagrange multipliers (Afe). 

In other words, the 5 kinematic degrees of freedom are made consistent with their 3D counterparts. The plate 
part of the 3D displacement, like warping in Section 2.3, is defined using averages weighted by (r^ • n). The Lagrange 
multipliers (A*,) are the magnitudes of the normal forces and moments due to the reaction of the 3D domain to the 
prescribed displacements: 


(17) 


fcG|l,5], f g{s,z) ■ n{s) ■^f{z)dz = Xk{s) f z'^iz) ■ n{s) ■ ^f.{z)dz = Xk{s) 

J Z J Z 


(18) 


The information concerning the last three stress components | 6 , 8 ]|, which involve only tangential action, is given by 
a classical Neumann condition. 

This Lagrangian coupling has the advantage of being applicable to zones of interest with no Dirichlet condition. 
It also makes the plate’s stress components in the 3D zone of interest directly available, which is useful in certain 
conhgurations of the global correction step (typically when no buffer zone is used). Unfortunately, this coupling 
requires the local orientation of the interface to be taken into account, which makes it hard to implement for curved 
domains. (In the above equations, (j) depends implicitly on s.) Also, it is difficult to choose a proper space for multiplier 
A( 7 c) in the case of nonconforming meshes: the LBB condition must be satished (see Section 3.5.2). 


3.1.3 Pure displacement descent 

Pure displacement descent is a very easy coupling to set up. The principle is to transmit as much information 
as possible using only Dirichlet conditions at the interface. This eliminates the potential instability problems of 
Lagrangian coupling and makes the coupling independent of the local orientation of the interface. The prescribed 3D 
displacement at the interface is given by: 

u{s,z) = ^ ®fe( 2 )Vfe(s) ^ Wj(z)Fj(s), s G 7 C, 2 : G [-h,h] (19) 

fc=l:5 i=l:8 

In other words, the classical plate displacements are enriched by warping vectors weighted by the generalized force 
components of the associated plate. 


3.2 Evaluation of the residual 

Now let us assume that a 3D problem has been solved through a Lagrangian or pure displacement descent step. This 
leaves us with a plate solution in w = w/ U ujb U loq and a 3D solution in U Db with kinematic continuity along 
7 c (in a sense which depends on the descent chosen). The residual (L^) is defined as the imbalance (in a plate sense) 
between the plate and the 3D domains along line 7 /: 

■ generalized plate forces along 77 viewed as the border of ujc U ujb 

]pb. 3 D(^) _ 3D generalized forces along 7 / viewed as the border of uji ( 20 ) 

Lfe = along 77 , fc G |1, 5] 

The residual is evaluated in a plate sense because the intention is to reuse it as a corrective load applied to the plate 
model. A 3D residual could be defined through a recovery of the plate forces, but, besides being unsuitable for the 
global/local strategy, it would be difficult to evaluate because an error of order (A^ is expected in the cr^z and ayz 

components and an error of order ^ 77 ^ is expected in the azz component. The fact that the 3D quantities have the 
correct orders of magnitude is interesting information which contributes to the validation of the hybrid model. 
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3.3 The global correction step 

The global correction step consists in solving the global plate problem after having updated the loading by residual 
(Rk) along interface 77 . The updating takes the following form: 



( 21 ) 


where Vf’* are the degrees of freedom of the plate test displacement field in the coordinate system of 77 . 

The iterations converge as long as the plate model of the zone of interest loi is stiffer than the 3D model D 7 
(classical fixed point). Otherwise, relaxation is necessary. Such an algorithm also converges for nonlinear monotonic 
problems [23]. 

In the linear case, the convergence can be improved, as is classically done in Schwarz methods, by using Krylov 
acceleration [50]. In this case, the full plate model (w) can be interpreted as a preconditioner to the hybrid model 
described in the following section. 

In cases where the zooming step has sufficient Dirichlet conditions {i.e. the Lagrangian and pure displacement 
cases), the bufferless limit case 77 = 7 c is acceptable. In the Neumann case, a nonzero buffer is necessary. In any 
case, overlapping is interesting with all methods because it reduces the small edge effects generated by the interface 
loading of the zooming step. 


3.4 Characterization of the limits 


At the convergence of the iterative process, the residual is zero and the hybrid model can be interpreted as follows: 

• In the zone of interest D 7 , it is a pure 3D model, 

• In the complementary zone wc, it is a pure plate model, 

• Along interface 77 , the stresses are in equilibrium, 

• Along interface 7 c, the plate’s kinematic continuity is satisfied (along with some other relations which depend 
on the zooming technique chosen), 

• In the buffer zone, the 3D problem in Ds is equivalent to the plate problem in lob in the following sense: 



:(St®e^))-Ff^) Nldx = Q 


Vfce [1,51, V<(>, 'iYl 


® n)) - = 0 on 77 , V/c G [1,5] '' ’ 

{u ■ {zf. ■ n)) - = 0 on 7 c, Vfc G [1, 5] 

In other words, the plate model obtained by integrating the 3D model through the thickness in Ds must be 
identical to the plate model in ujb- 

One should note that in all cases 3D displacements orthogonal to the Saint-Venant tractions are likely to develop in 
Db, but these should be evanescent and should have negligible effect on D 7 . 


3.5 Implementation into a finite element code 

3.5.1 Preprocessing of the Saint-Venant stresses and warping vectors 


This preprocessor requires an auxiliary 3D problem to be built with the same stacking sequence as the original problem, 
but a simpler geometry. Minimal Dirichlet conditions are prescribed, except for shear-dominated experiments (j = 5 
and j = 8 ). 

An important property is that because the plate’s kinematics (spanned by functions (Sfe)) is linear through the 
thickness, it is represented exactly by the 3D FE approximation. This means that the calculation of the plate’s 
generalized stresses can be performed directly on the nodal reactions (which are classical outputs of FE software). Let 
(F*) be the nodal reactions through the thickness calculated on a face with normal e^. If z* is the height of the 
node), one has: 


Nxx = 

II 

II 

II 

= 

II 

M' 

^xy — {zO'xy) - 



i 


i 


(23) 


Qx = {<Jxz) = ^FI 
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3.5.2 Issues related to the Lagrangian approach 

The implementation of the Lagrangian approach requires great care in choosing the finite element approximation 
spaces. 

In the case of matching meshes and interpolations (typically 8-node hexahedra on the 3D side and linear plate 
elements on the plate side), each plate node faces a column of 3D nodes across the thickness. A simple choice consists 
in using five multiple-point constraints (MFCs) per plate node and corresponding 3D nodes. 

The case of nonmatching meshes or interpolations raises the classical issues of mixed approaches: one must check 
the LBB-stability of the scheme with a well-chosen discretization of the Lagrange multiplier field along the interface. 

For example, for a linear plate element and a bilinear 3D element (20-node hexahedron), one method is to use two 
plate elements facing one column of 3D elements (so each plate node faces a column of 3D nodes), then use a constant 
Lagrange multiplier for each column of 3D elements. 

4 Numerical study of the various zooming steps 

In order to illustrate the descent step, let us consider the cantilever plate of Figure 12, with L = a = 5 mm, h = 0.2 
mm and thickness 2h. The plate consists of four orthotropic plies (—45/45)^ with the following mechanical properties: 

El = 25 MPa 
El = Epf = 1 MPa 
Glt = Gln = 0.5 MPa 
Gtn = 0.2 MPa 
vlt = vtn = vln = 0.25 



The prescribed displacement at the edge of the plate is Ud • = 2h. The zone of interest is the clamped side 

because this is where 3D localization takes place. The 3D model was meshed with quadratic hexahedral HEXA20 
elements (20 nodes). A regular mesh was used with 4 elements per ply in the z-direction and 20 elements in the x 
and y directions. The global model defined over the entire plate was obtained using Dhatt and Batoz homogenization 
[3] and a discrete shear formulation. The reference is a full 3D model over the whole domain. The calculations were 
carried out using Python scripts interfaced with Code_Aster. 

Since the zone of interest is clamped, it was possible to use all three zooming techniques. Figure 13 shows the 
evolution of some components of the 3D stress along the x axis. The reference is defined in D (a; G [0,10]) while the 
3D parts of the hybrid models are defined in zone x G [0,5]. In all three methods, as expected theoretically [1], axz 
has approximately the same order of magnitude as and azz has approximately the same order of magnitude 

as (Jxxi^)- One can observe that in all cases a perturbation appears near the interface (a: = 5). However, since the 
generalized forces transmitted at the interfaces are correct, this local effect decreases rapidly and vanishes within a 
2h band. The 3D solutions are close to the reference solution. Altogether, the perturbations caused by the warping 
technique seem smaller than those caused by the Neumann technique, which, in turn, are smaller than with Lagrangian 
coupling. The presence of these undesirable stress concentrations is what triggered the extension of the zone of interest 
by a 3D buffer zone in order to dampen the artificial edge effects. 

Figure 14 shows the evolution of the generalized stresses along the interface for the initial plate calculation and 
for the 3D calculations. By construction, Neumann zooming produces the same generalized stresses as the plate 
calculation. Conversely, the other zooming techniques lead to slightly different results, which justifies the iterations 
between the local and global models. As expected, the differences between the plate model and the zooms are 
particularly significant near the edges. 
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(b) Normal stress 



(c) Peeling stress (d) Shear stress 

Figure 13: Evolution of stresses along the x-axis 

5 Evaluation of the hybrid model based on warping 

In this section, we study the application of the iterative coupling technique in order to achieve convergence toward 
the hybrid model. Based on previous experiments, we chose the warping-based descent algorithm (see 3.1.3) and used 
a nonzero buffer zone (see 3.3). 

The objective was to model the holed composite plate shown in Figure 15. The local model consisted in a 3D 
representation of the hole and its vicinity with enough extent to capture the 3D edge effect. The global model was a 
plate, in which we could choose to represent (or not) the hole in the zone w/ corresponding to the 3D model. We used 
the same orthotropic layers as in the previous section. The dimensions were L = a = 60 mm and h = 1 mm. The 
radius of the hole was r = 2 mm. The loading was a transverse displacement Ud = 0.45 along the right side. 

The zone of interest w/ was a square of side 20 mm and the width of the buffer zone ojb was 2 mm. It was meshed 
with 24 elements per side and 16 elements in the thickness. 

5.1 Convergence toward the hybrid model 

Let us define the relative norm rji of the residual L of Section 3.2 as follows: 

2 _ f-yi Ar Ar ^ 

Vi z~p T ~ p T ~ p 

max-yj • np max.^, |• nP max.^, \Q^ • nP 

Figure 16 shows the decrease in rji as a function of the number of iterations for the classical fixed point algorithm and 
for the conjugate gradient algorithm. The complexity of the plate model in the zone of interest w/ is a parameter of 
the method which influences the convergence, but not the limit. We considered two cases, depending on whether the 
hole was represented in the plate model luj or not. 

When the hole was included, the plate model gave a better representation of the zone of interest, leading to much 
less initial imbalance. In both cases, the convergence was fast. In particular, one can observe that two conjugate 
gradient iterations of the model without a hole led to a similar quality approximation as the classical submodeling 
technique applied to the holed model (which is more difficult to generate and to mesh). 
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(b) Bending moment 




y y 


(c) Torsion moment 


(d) Normal shear 


Figure 14: The generalized forces along the interface 



Figure 15: The coupling of the two models (left) defines the hybrid model (right) 



Iteration 


Figure 16: Evolution of the relative norm of the residual as a function of the number of iterations 


5.2 Validation of the hybrid model 

Now let us study the hybrid model obtained after a sufficient number of iterations of the coupling procedure and 
compare the results of this model with those of a full 3D calculation of fl in terms of the reliability of the 3D quantities 
obtained. 
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The quantities of interest chosen are the stresses in the vicinity of the hole, which would initiate delamination. 
Figures 17(a,b,e) show the variations of some stress components along the particular line described in Figure 17f. A 
comparison of the full 3D reference values, the values using submodeling (which is equivalent to the hybrid model 
with a single descent and no iteration) and the values given by the hybrid model after 5 iterations clearly shows that 
submodeling is very unreliable while the hybrid model is almost identical to the reference model. 

Figures 17(c,d) show the evolution of the generalized forces along the same line. A comparison of the reference 3D 
calculation, the hybrid model at iterations 0 and 5 and a direct (holed) plate model shows that both the plate and 
the submodel lead to poor results while the hybrid model comes very close to the reference. 

Table 2 quantifies the changes induced by the iterations in the plate quantities in the plate domain: the final values 
of the deflection, shear stress and moment {i.e. the values corresponding to the hybrid model) at a point of 7 / (x = 20 
mm, y = 0 mm) are compared with the initial values obtained with the plate model (with or without a hole). Of 
course, the model without a hole had significant errors, but even the model with a hole was improved by taking into 
account 3D phenomena. (In particular, the shear stress was underestimated by more than 5%.) 


— 


with a hole 

without a hole 

error in w 

(%) 

-1.93-10-2 

-1.97 

error in 

(%) 

5.99 

-35.94 

error in M^x {%) 

1.18 

3.3 


Table 2: The corrections between the first descent and the fifth iteration 


6 Conclusion 

In this paper, we proposed a new strategy for modeling thin composite structures by coupling a global plate model 
with a refined 3D model in zones of interest. The objective is to combine the efficiency of plate modeling for most of 
the structure with the added accuracy of 3D calculations where necessary (near edges, holes or defects). The hybrid 
plate/3D model is obtained after the convergence of a nonintrusive iterative process which enables one to use standard 
commercial software and industrial plate models without even having to represent the critical zones. The coupling 
is achieved by means of a preliminary analysis in which numerical Saint-Venant stress and warping fields which are 
suitable for the stacking sequence and the chosen discretization are recovered automatically. Several coupling methods 
are possible, but the warping-based technique is the simplest to implement and seems to lead to the least amount 
of perturbation at the interface. The convergence can be monitored by looking at the imbalance (in a plate sense) 
between the 3D domain and the plate domain. In the case of linear 3D domains, it can be accelerated by using a 
Krylov solver. 

Our works in progress concern the application of the method to the simulation of complex assemblies of composite 
plates. In this case, the global model uses classical simplified ID connectors to maintain computational efficiency while 
realistic 3D nonlinear models of bolts, rivets and damageable composites are used in the local models. A particular 
issue arises from the use of different in-plane discretizations at the interface between the local model and the global 
model. Two-scale approaches, such as that of [26] in the context of the XFEM, have interesting advantages, such as 
the possibility to transfer average static plate quantities to the 3D model without introducing local forces and spurious 
edge effects. 
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